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Abstract. The two-dimensional layout optimization problem reinforced by the efficient space 
utilization demand has a wide spectrum of practical applications. Formulating the problem as a 
nonlinear minimization problem under planar equality and/or inequality density constraints, we 
present a linear time multigrid algorithm for solving correction to this problem. The method is 
demonstrated on various graph drawing (visualization) instances. 
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1. Introduction. The optimization problem addressed in this paper is to find 
an optimal layout of a set of two-dimensional objects such that (a) the total length 
of the given connections between these objects will be minimal, (b) the overlapping 
between objects will be as little as possible, and, (c) the two-dimensional space will 
be well used. This class of problems can be modelled by a graph in which every vertex 
has a predefined shape and area and each edge has a predefined weight. While the 
first two conditions are straightforward, the third requirement can be made concrete 
in different ways. To see its usefulness, consider for example, the problem of drawing 
the "snake" -like graph shown in Figure [T]^ a). Most graph drawing algorithms would 
draw it as a line or a chord. In that case, when the number of nodes is big, the space is 
used very inefficiently, and the size of the nodes must decrease. One possible efficient 
space utilization for the graph " snake" is presented in Figure [l|b) . 




(b) 



Fig. f.f. Possible ways to draw the "snake"-like graph: (a) when the drawing area is not used, 
the size of the nodes must decrease; and (b) a clearer picture is obtained when the space is used 
efficiently. 

In many theoretical and industrial fields, this class of problems is often addressed 
and actually poses a computational bottleneck. In this work we present a multilevel 
solver for a model that describes the core part of those applications, namely, the 
problem of minimizing a quadratic energy functional under planar constraints that 
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bound the allowed amount of material (total areas of objects) in various subdomains 
of the entire domain under consideration. 

Given an initial arrangement, the main contribution of this work is to enable a 
fast rearrangement of the entities under consideration into a more evenly distributed 
state over the entire defined domain. This process is done by introducing a sequence 
of finer and finer grids over the domain and demanding at each scale equidensity, 
that is, meeting equality or inequality constraints at each grid square, stating how 
much material it may (at most) contain. Since many variables are involved and since 
the needed updates may be large, we introduce a new set of displacement variables 
attached to the introduced grid points, which enables collective moves of many orig- 
inal variables at a time, at different scales including large displacements. The use of 
such multiscale moves has two main purposes: to enable processing in various scales 
and to efficiently solve the (large) system of equations of energy minimization under 
equidensity demands. The system of equations of the finer scales, when many un- 
knowns are involved, is solved by a combination of well-known multigrid techniques 
(see 131 SI I14|). namely, the Correction Scheme for the energy minimization part 
and the Full Approximation Scheme for the inequality equidensity constraints defined 
over the grid's squares. We assume here that the minimization energy functional has 
a quadratic form, but other functionals can be used via quadratization. The entire 
algorithm solves the nonlinear minimization problem by applying successive steps of 
corrections, each using a linearized system of equations. 

Clearly, for each specific application, one has to tune the general algorithm to 
respect the particular task at hand. We have chosen here to demonstrate the perfor- 
mance of our solver on some instances of the graph visualization problem showing the 
efficient use of the given domain. Let us review a few applications that have motivated 
our research. 

Graph visualization addresses the problem of constructing a geometric repre- 
sentation of graphs and has important applications to many technologies. There are 
many different demands for graph visualization problems, such as draw a graph with 
a minimum number of edge crossings, or a minimum total edge length, or a predefined 
angular resolution (for a complete survey, see [2J). The ability to achieve a compact 
picture (without overlapping) is of great importance, since area-efficient drawings are 
essential in practical visualization applications where screen space is one of the most 
valuable commodities. One of the most popular strategies that does address these 
questions is the force directed method [7 which has a quadratic running time if all 
pairwise vertex forces are taken into account. There are several successful multilevel 
algorithms [10^ developed to improve the method's complexity. However, reducing 
the running time in these models usually means a loss of information regarding those 
forces. 

Representation of higraphs. Higraphs, a combination and extension of graphs 
and Euler/Venn diagrams, were defined by Harel in [9j. Higraphs extend the basic 
structure of graphs and hypergraphs to allow vertices to describe inclusion relation- 
ships. Adjacency of such vertices is used to denote set-theoretic Cartesian products. 
Higraphs have been shown to be useful for the expression of many different semantics 
and underlie many visual languages, such as statecharts and object model diagrams. 
The well-known force-directed method has been extended to enable handling the vi- 
sualization of higraphs [8^. For small higraphs it has indeed yielded nice results; but 
because of its high complexity, it poses efficiency challenges when used for larger 
higraphs. 
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Facility location problem. In this class of problems the goal is to locate a 
number of facilities within a minimized distance from the clients. In many industrial 
versions of the problem there exist additional demands such as the minimization of the 
routing between the facilities and various space constraints (e.g., the factory planning 
problem) while given a total area on which the facilities and clients could be located 
(for a complete survey, see [6 ). 

Wireless networks and coverage problems have a broad range of applications 
in the military, surveillance, environment monitoring, and healthcare fields. In these 
problems, having a limited number of resources (like antenna or sensor), one has 
to cover the area on which many demand points are distributed and have to be 
serviced. In many practical applications there are predefined connections between 
these resources that can be modeled as a graph [T^l |5J [TT] . 

The placement problem. The electronics industry has achieved a phenomenal 
growth over the past two decades, mainly due to the rapid advances in integration 
technologies and large-scale systems design - in short, due to the advent of VLSI. 
The number of applications of integrated circuits in high-performance computing, 
telecommunications, and consumer electronics has been rising steadily, and at a very 
fast pace. Typically, the required computational power of these applications is the 
driving force for the fast development of this field. The global placement is one 
of the most challenging problems during VLSI layout synthesis. In this application 
the modules must be placed in such a way that the chip can be processed at the 
detailed placement stage and then routed efficiently under many different constraints. 
This should be accomplished in a reasonable computation time even for circuits with 
millions of modules since it is one of the bottlenecks of the design process. For a most 
recent survey of the placement techniques see PTS^. 

The paper is organized as follows. The problem definition is described in Section 
|2] The multilevel formulation and solver are presented in Section |3] Examples of 
graph drawing layout corrections are demonstrated in Section [4] 

2. Problem definition. Given a weighted undirected graph G — {V,E), let 
v{i) > be the (rectangular) area of vertex (node) i G V, i = 1, \V\, and Wij the 
non-negative weight of the edge ij between nodes i and j {wij = if ij ^ E) . Also, 
assume a twodimensional initial \a.yo\it is given; that is, the center of mass of node i is 
considered to be located at (xi,yi) within a given rectangular domain. The purpose 
of the optimization problem we consider is to modify the initial assignment {x^y) by 
{^XT^y) so as to minimize the quadratic functional 

'^{5^x.5y) w^j{{x, + 5^^^ - ij - 5^^f + {y^ + 5y^ ~ Vj ~ 5y^f) , (2.1) 

subject to some equidensity demands on the area distribution of the nodes within the 
given rectangle. To apply such constraints, we discretize the domain by a standard 
grid Q consisting of a set of squares S{G), where each square s € S{G) is of area 
A — hxhy and and hy are the mesh sizes of G in the x and y directions, respectively 



(see Figure 2.11. Denote by T(s) the total area of the vertices overlapping with the 



square s; that is, T(s) is the sum over all the nodes coinsiding with s, each contributing 



the (possibly partial) area that overlaps with s (see Figure 2.2). 

The planar constraints (i.e., the constraints that are distributed over the 2D 
plane, where each constraint defines a demand regarding some bounded area) can 
then simply state how much area is required to be in every square; that is, for each 
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Fig. 2.1. Example of a grid Q with 25 grid points and 16 squares. The grid points and squares 
are labeled by pi and bold numbers, respectively. 
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Fig. 2.2. Example o/T(s) for square 6. The total area of vertices overlapping with square 6 is 
dashed. 



square s € S{Q) the constraint is either T(s) = M{s), or T(s) < M{s), where M(s) 
is the amount of nodes area desired or allowed for square s. 

The constrained optimization problem with equality or inequality formulation can 
thus be summarized by the following 



minimize £ (given by (2.1 1 



subject to T{s) = {<)M(s), ys £ S{g) 



(2.2) 



3. The multilevel formulation and solver. The aim of the current work is to 
provide a fast first-order correction to the given approximate solution; that is, we are 
looking for such a displacement that would in some optimal sense (to be defined below) 
improve the planar equidensity demands and/or decrease €. (Note that unconstrained 
minimization of 2; will bring all nodes to overlap at a single point, and thus we may 
often observe an increase in (S upon removing some of the initial overlap.) 

To enable a direct use of the multigrid paradigm, and motivated by the need to 
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perform collective moves of nodes (as explained in the introduction) , we have actually 
reformulated the problem (2.11 as described in Section 3.1 The multilevel solver of 



the (reformulated) system (3.101 below is introduced in Section 3.2 This system of 



equations actually has to be solved for a sequence of different grid sizes to enhance 



the overall equidensity for a variety of scales as presented in Section 3.3 



3.1. Formulation of the correction problem. We have first introduced two 
new sets of variables u and v that correspond to displacements in the horizontal and 
vertical directions, respectively. These variables are located at the grid points V{Q) 
which are sequentially counted from to |7'(t/) | — 1 as shown in Figure 2.1 Each point 



P G T-'iQ) is associated with two variables Up and Vp that influence the displacements 
of all the nodes located in the (up to four) squares intersecting at p. For example, 
the horizontal update of (the center of mass of) node j, depicted in Figure 2.1 is 
obtained from points pi2,Pi3,Pn and pig: 



Xj ^ Xj + ai2jMi2 + ai3jWi3 + anjuij + aisjWis, 

where ai2, ais, (Xn and aig, are the standard bilinear interpolation coefficients (their 
sum equals 1). The vertical coordinate yj is updated from the v variables using the 
same coefficients. 

For a node i denote the set of four closest points in 1^(0) (the corners of the 
square its center of mass is located at) by c(i). The new quadratic energy functional 
we would like to minimize for u and v given a current layout {x,y) of G (i.e., the 
coordinates of node i are initialized with (xi,yi)) is 



(£(m, v) 



ijeE ^ pec(j) P<^c(j) p€c{i) P&c{j) 



E 



(3 



where api are the bilinear interpolation coefficients. 

The reformulation of the equidensity constraint in terms of the displacement 
variables relies on the rule of conservation of areas. The initial total amount of vertex 
areas at each square equals the current actual amount of areas dictated by (i,y). To 
estimate the amount of areas flowing inside and outside a given square induced by 
the u and v displacements, we assume the nodes are evenly distributed inside the 
squares. Under this assumption it is easier to estimate the amount of area being 
transferred between two adjacent squares as explained below. Consider, for example, 
a square s. Denote by T^^; ^ ;,-)(s) the total area of nodes overlapping with its right 
(left, top, bottom) neighbor square. Let Uj.t{rb,it,ib){s) be the u values at the right-top 
(right-bottom, left-top, left-bottom) corner of s as shown in Figure 3.1 To estimate 



the amount of areas entering s from the right we flrst calculate the average area (per 
squared unit) in both squares: (T(s) + Tr{s)) /2A. We have to multiply this by the 
actual entering area (of nodes), which is a rectangle of height hy, the length of the 
border between the two squares, and width, which is the average of the u displacement 
at the middle of that border, namely, {urt + Urb)/2. Thus the overall contribution of 
area from the right is approximated by 



2A 



■ K, 



Urtis) + Urbis) 
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A similar term is calculated at the left, and with v instead of u also at the top and 
bottom. Note that if the assumed direction of flow is wrong the resulting displacement 
will just turn out to be negative. 

The entire constraint for a square s stating that the net flow of areas into the 
square should be equal to or be smaller than some demand M(s) minus the current 
area in u, is given below: 

T{s)+Tr{s) Urt(s) + Urb(s) T(s)+T;(s) Uit (s) + Ulb («) , 

'"^"^'^ ^ 2A 2 2A 2 + 

T(s) + Tt(s) i^rt(s) + i'it(g) T(s) + Tb(s) l'rb(s) + vihjs) , . , . 

2A 2 2A 2 - ^'^^'^ ~ ■ 

(3.2) 

Next, to enforce the natural boundary conditions on u and v, namely, to forbid flows 
across the external boundaries, we simply nullify all corresponding Up on the right and 
left boundary points Bu{Q), and Vp on the bottom and top boundary points Bi,{Q). 
Then the entire constrained optimization problem in terms of u and v and the initial 
approximation (x, y) is given by 



minimize u) (given by (3.1 ) 

subject to eqO(s) = (<)M(s) -Tfs) , Vs e S{g) ; 

li p Bu{g) i'hcn Up = Q ; 
if p e Bv{g) then Vp = . 



We will simplify the formulation of ( |3.3| by the concatenation of the two vectors u 

and V into one u = [{ui}-^Q | {wi}i=o ]■ ^^^^ omit the boundary 
conditions by directly replacing all variables in Bu{G) ^ Bv{g) by 0, and so, from now 
on, we will refer to <B as 

&{u) ^^^q^jUiUj + ^kUi + C , (3.4) 

where i,j run over all the indices in u\ BuiG) \ By{g), C is a constant and q^, li are 
the coefficients calculated directly from the previous deflnition (3.11 of £. Similarly 
rewrite each eqc)(s) in (3.2) as 

eqO(s) = ^ QsiUi = {<)bs, (3.5) 

i 

where bg — M{s) — T(s). 

Denote by As, s € S{g) the Lagrange multiplier corresponding to the equiden- 
sity constraint of square s. If all the constraints are equality ones, the Lagrangian 
minimization functional is 

£(u,A) = £(u)+ J2 As(eq^(s)-fes) • (3.6) 

ses{g) 

So, we are looking for a critical point of the Lagrangian function, which is expressed 
by the system of linear equations 



V£(u,A) 



Vu£(u,A) 
Va£(u,A) 



. (3.7) 
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There are at least two factors that may cause (3.7 1 to be singular. First, the rank 



of Vi2(u, A) is always less than its size by at least 1. This arises from the equations 



of equidensity constraints in (3.7): their sum always equals zero. The reason is that 



under the boundary constraints the total amount of in-flows is always equal to the 



total amount of out-flows. In fact, the second summand in (3.6) can be replaced by 
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Fig. 3.1. The horizontal direction flows of area considered for the square s (colored by gray) in 
the equidensity constraint 



J2 iXs + Z){cqHs)-bs) 

sesiG) 

for any Z without changing the minimization of £ since 

sesiG) 

Thus, important are not the values of As but only their differences, and the singularity 
can be treated by an additional constraint, say, dsXs = 0, where = 1 Vs G 
(the introduction of dg is necessary for the recursion of the multilevel solver; see 
Section 3.2.1). The additional term in £(u, A) is rjJ^s^sXs, where is a "pseudo- 



Lagrange" multiplier. The following proposition (with fc = 1) motivates the non- 
singularity of £ with J2s d-sXs ~ 0. 

Proposition 3.1. Given a symmetric nxn matrix A, for which rank(A) — n—k, 
let Xi, « = 1, he an orthgonal basis oj the null space of A, that is, Axi = 0. Then 
the following block matrix B is nonsingular 



B 



A 


X 








where X — {xi, X/t) is an n x k matrix of rank k 



Proof. Let y be any vector in M"+'''. Denote by y' the first n components of y 
and by y" the last k components, that is, y = (f")- ^^^^ prove that if By — 0, 
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then y — 0. The vector By can be written in the following block form: 

Ay' + Xy" 



By = 



Multiplying Ay' + Xy" = by from the left implies that y" — 0, and hence 



Ay' — and y' ~ Q^i^i- Substituting the last relation into each of the last k 



rows of B implies x^y' — xj J2i=i '^i^i — '^j^J^j — thus aj = ioi j — 1, fc 
yielding y' = 0. Since y' — and y" = we may conclude that y = as needed. □ 



The second kind of singularity in (3.7 1 may appear from possible empty squares. 



This can be treated by adding a summand to (3.6 1 that minimizes the total sum of all 
corrections (iJ2i that is, adds a 2/3-term to the diagonal of Vu£, where /3 is small 
enough to cause only negligible change in a solution. This will prevent the inclusion 
of zero-rows in Vu£, while possibly also bounding the size of each correction in the 
solver below. 

To summarize, the pseudo-Lagrangian functional £ for our correction problem 
with equality constraint is 



(3.8) 



2 

■ij » i ses{g) i sG5(e) 



leading to the following system of equations 

J2i{aszUi - bs) + r/ds^O , Vs G 5(0) 

J2seS{g) ^s^s = . 

(3.9) 

Since in real world situations the total area is usually bigger than the total area 
of all the vertices, the redefined minimization problem under inequality constraints 
will generally have the form 

minimize 2;(u) (given by (|3.4 l)) , , 

subject to eqc)(s) < bg , Vs G S{Q) (given by (|3.5f ) . \ ■ ) 



3.2. Multilevel solver for problem ( 3.10[ ). To solve the constrained mini 



mization problem (3.10), we use multigrid techniques: standard geometric coarsen- 
ing, linear interpolation. Correction Scheme for the energy minimization and the Full 
Approximation Scheme for the equidensity inequality constraints; all are presented in 
Section [3.2. 1[ In addition, we have developed a fast window minimization relaxation 



as explained in Section 3.2.2 The multilevel cycle is schematically summarized in 



Section 3.2.3 in Algorithm 2D-layout-correction. 



3.2.1. Coarsening scheme. When the geometry of the problem is known we 
can choose a coarser grid by the usual elimination of every other line, as shown in 
Figure |3.2[ The correction computed at the coarse grid points will be interpolated 
and added to the fine grid current approximation. Let us introduce the notation 
distinguishing between fine and coarse level variables and functions. By lowercase and 
uppercase letters we will refer to the variables, indexes, and coefficients of the fine 
(ui, i, Qj, etc.) and the coarse (U/, /, Qj, etc.) levels, respectively. The subscripts / 
and c will be used to describe the energy (Bf and €c and pseudo-Lagrangian £/ and 
£c functions at the fine and the coarse levels, respectively. 
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Fig. 3.2. Geometric coarsening. The equidensity constraints of every four similarly patterned 
squares at the fine level form one equidensity constraint at the coarse level. 



Thus, the minimization part of the pseudo-Lagrangian (3.8 1 at the fine level is 

+ V/zU, . (3.11) 



-E 

^j 



(Note that we have omitted the /3 term from the following derivation since it is merely 
an artificial added term.) Given a current approximation u of the fine level solution 
u and a correction function U calculated at the coarse level variables U, u will be 
corrected by 



U; ^ Ui + ^ fti/U/ 



(3.12) 



I3i 



where the notation means that the sum is running over all coarse gridpoints pi 

from which standard bilinear interpolation is made to the fine gridpoint pi. 

Expressing the fine level energy functional Ef in terms of the coarse variables by 
substituting ( |3.12[ ) into ( |3.1ip yields 



ij I3i 



J3J 



E 

I3i 



a^/U/) 



1 



2 ^ 

I, J 



QjjlJilJj + Y,LiUi + C , 



where Qu = J2iei QijC^uajj, Lj j Qiji»-j<^ii + J2iei ^i^n C is a constant. 

Thus, the coarse level energy functional will be of the same structure as the fine level 
one, namely. 



IJ 



For each fine square s the equidensity constraint eqfi(s) is given by (3.51. The 
coarse equidensity constraints are constructed by merging 2x2 fine squares into one 
coarse square S. The expression "s G S"' will refer to running over the four fine 



squares s that form the coarse square S (see Figure 3.2 1. The S'-th planar equidensity 
constraint of the coarse level (in the case of equality constraints only) is obtained 
again by the substitution of ( [3.12^ : 



ses i 
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where Asi = Y^ic^i J2 



.sc^s^sictii and Bs = Z]s6s(^s ~ Y.iO-siVLi). Similarly (in the 
case of equality constraints) , the additional ry-constraint over all squares at the coarse 
level as inherited from the fine level is DsAs = 0, where Ds = J2ses 

To complete the description of the coarse equations, we still need to transfer the 
equidensity inequality constraints. For this purpose we will use the Full Approximation 
Scheme (FAS), which is the general multigrid strategy applied to nonlinear problems 
(see [31 m [T3]). In fact, there is no need for the FAS for the equality equidensity 
constraints since it is a linear problem that can be solved by the regular Correction 
Scheme (CS). The FAS-like coarsening rules are needed and applied only on the 
set of equations derived from the equidensity inequalities. Thus, our scheme is a 



combination of the correction scheme for the energy equations derived from (3.111 



and (3.121 and FAS-like rules for the equidensity equations. 



To derive these equations we need to calculate the residuals for both the fine and 
coarse grids. If £/ is the pseudo-Lagrangian of the fine level system defined by 



s) +??y^rfsAs 



(3.13) 



where is given by (3.11), then the u^-th residual of V£/, where As is the current 



value of the Lagrange multiplier As, is 



Thus, the residual corresponding to the variable XJj of V£c (where V£c is the coarse 
level system of equations analogous to ( |3l3| ) IS 



(3.14) 



where an are as in (3.121; that is, the fine-to-coarse transfer is the adjoint of our 



coarse-to-fine interpolation. The residual of the s-th equidensity constraint is 

rl'^^ = bs-^ asiUi - fids , 



where s runs over all fine squares and ry is the current value of 77. Therefore, the coarse 
equidensity residual of square S is 



(3.15) 



Finally the residual of the ry-constraint is 

r,, = -^dsAs = Rh- 

s 

Denote by LP{I) the linear part of the Uj-th equation in the system V£c 
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From the FAS rule for the I-th coarse equation stating that LP{I) ~ Rf- 
approximation of LP{I) , we can derive the /-th V£c equation 



the current 



0, 



(3.16) 



where Rj is given by (3.141. U' 



and 



\ Esss ^s- Similarly, the S-th. square 



coarse equation for the equality (inequality) constraint is 



AsiVi + HDs 







H^'Ds = (<)0 , 



where R'^g^ is given by (3.15). The last equation for the _ff-constraint is 



E DsKs -Rh-J2 



0. 



(3.17) 



(3.18) 



Note that equations (3.161 to (3.18) are the coarse grid equations analog to the system 



(3.9). (A 2/3U/ term may be added to (3.16) for stabihty if needed.) The correction 



received from the coarse level for the u variables is given by ( 3.12 ) and for the Lagrange 
multipliers A by 



A. 



Xs + A 



S3s 



(3.19) 



3.2.2. Relaxation. In our multigrid solver, as usual, the relaxation process is 
employed as the smoother of the error of the approximation, before the construction of 
the coarse level system and immediately after interpolation from the coarse level. For 
this purpose we have developed the Window relaxation procedure, which extracts from 
the entire system small subproblems of m x m squares and solves each separately, as 
explained below. The running time of the entire relaxation process strongly depends 
on the algorithm for solving one window. There exist many versions of well-known 
algorithms for the quadratic minimization problem under linear inequality constraints 
(for a survey see [1 ). However, since each window need be solved only to a first 
approximation (because of the iterative nature of the overall algorithm) , in order to 
keep the running time low, we have implemented a simple algorithm for approximately 
solving each single window, as presented in SingleWindowSolver. 

Let W = {s G iS(fJ)| all squares within an m x m super-square} be a window of 
squares. To solve the quadratic minimization problem in W, we fix at their current 
position all u outside W, as well as all those that are on the boundary of W and 
represent movement perpendicular to the boundary. The minimization is done under 
the set of equidensity constraints for the squares s € W. The solution process for each 
single window is a simplified version of the active set method and is iterative. At each 
iteration t, for a given u we first extract the set (denoted by St) of squares for which 
the respective inequality equidensity constraints are violated or almost violated: 

St = {seW\ zqd{s) >b,-e} , 

where e is positive and sufficiently small but not too small to make St numerically 
unstable (we have used e = 0.0001* (the square's area)). Then the inequality con- 
straints of St are set to equalities ignoring the other inequality constraints. Let Vw 
be the set of all displacement indexes inside W (including those on the boundary of 
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yV directing parallel to it). For every u^, i G Vw we associate a correction variable 6i 
and we reformulate the pseudo-Lagrangian for W as a functional of the Si variables 
as follows: 

where is the current value of and the f3 term is added for stability with P — I. 
Solving V£vy((5, A) = we obtain the corrections 6i for Ui,i £ Vw, which confine the 
respective active set variables to the boundary of the equality constraints manifold. 
However, while accepting this correction we may violate other inequality constraints 
that were already satisfied at the previous iteration i — 1. Let us call this set of 
new unsatisfied constraints St- One way to overcome this problem is to accept only 
a partial correction 9Si, i G Vw, where 9 is the smallest number that brings some 
constraint from St to equahty. Accepting the correction 9Si does not violate any 
constraint from St- At this point we accept this partial correction and continue to 
the next iteration t+1, excluding from the redefined St the set of satisfied (by equality) 
constraints from St with negative Lagrange multipliers A.,. 



SingleWindowSolver(yV, u) 
begin 

t = 

Repeat until " optimal enough" (explained at the end of Section H 
If i = 

St — {the violated equidensity constraints} 

Else 

St = {the violated equidensity constraints}\ 

{those from iteration t ~ 1 which satisfy equality and have A^ < 0} 
Solve V£vv(i5; A) = and extract the smallest 9 
Accept the correction u ^ u + 6'(5 

t^t+1 

end 



To achieve corrections for all variables, we will cover by these windows the entire 
area in red-black order |14j . For computational reasons we have chosen to apply this 
relaxation for very small windows (of size 4x4 squares) . To minimize the effects of 
the boundary constraints in the windows and to enforce the equidensity constraints 
over different super-squares, we scan the entire domain two more times: once with 
half-window size shift in the horizontal direction and once in the vertical. Thus the 
overall relaxation process covers the domain three times. 

3.2.3. The multilevel cycle. Having defined the window relaxation, the inter- 
polation, and the coarsening scheme, the multilevel cycle naturally follows. Starting 
from the given approximation (i, y), discretize the domain by a standard grid on 



which the u variables are initially defined. Construct the system of equations (3.9 1, 
and solve for the u variables as follows. After applying i>i window relaxation sweeps, 
define the coarser level equations for the coarser grid, apply i^i window relaxation 
sweeps there, and continue to a still coarser level. This process is recursively repeated 
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until a small enough problem is obtained. Solve this coarsest problem directly, and 
start the uncoarsening stage by interpolating the solution of the coarse level to the 
finer levels followed by V2 window relaxation sweeps on the finer level. Repeat until 
the correction to the original problem is obtained. This entire multilevel cycle, usually 
referred to as the V-cycle, is summarized in procedure V-cycle-correction below, 
where the superscript index refers to the level number. (We have used vi = V2 — 'i) ■ 



V-cycle-correction(^*, \i\ C\ X\ V£') 
begin 

If is a small enough grid 

Solve the problem exactly 
Else 

Set u' = 

Apply vi Window relaxation sweeps 
Construct the coarse level grid 

Define C*+^ to be the set of equidensity constraints 

Initialize the system of equations V£*"'"^ given by (3.16 1-(3.18 1 
Initialize u*+^ and A'"*"^ 

V-cycle-correction(^;*+i, u'+i, C'+i, A*+\ V£^i) 
Interpolate from level i + 1 to level i using (3.12) and (3.191 
Apply V2 Windows relaxation sweeps 
Return u' 



3.3. The Full MultiGrid external driving routine. The solution of (3.9 1 
is primarily dependent on the chosen grid size. To enforce equidensity at all scales, 
it can be used within the Full MultiGrid (FMG) framework. This is done by using 
a sequence of increasing grid sizes (progressively finer meshsizes), while employing a 
small number of V-cycles for each grid size. We emphasize that the original problem 
(2.2) is highly nonlinear, while the system of equations with corrections in term of 
the displacement u is linearized around the current solution {x,y). Therefore, only a 
small correction should actually be taken from the u displacements when these (i,y) 
are being updated. Then, a new linear system can be formulated around the new 
solution to obtain a new correction, and so forth. Thus, by small steps of corrections 
we solve the original nonlinear problem via the corrections calculated from the linear 
system of equidensity constraints. For instance, we have tried to employ grids of sizes: 
2, 4, 8, ... up to a grid with number of squares comparable to the number of nodes 
in the graph. For each grid size the corresponding set of equations (in terms of the 
displacement u) is solved either directly (for small enough grids) or by employing 
the V-cycles described in Section |3.2.3[ In either cases the obtained solution u is 
interpolated back to the {x,y) variables, introducing the desired correction to the 
original variables of the problem. Various driving routines can be actually used: each 
chosen grid size may be solved more than once (e.g., use grids 2, 2, 4, 4, 8,...); the 
entire sequence of grids may be repeated (e.g., 2, 2, 4, 4, 8,... , 2, 2, 4, 4, 8,...), and 
so on. (see Section |4] for examples). These parameters should in fact be optimized for 
each application according to the concrete needs of the model. 

The entire algorithm for the two-dimensional layout correction is summarized 
below in Algorithm 2D-layout-correction, where the superscript refers to the 
current chosen grid size. 
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2D-layout-correction(graph G, current layout {x,y)) 
begin 

Apply for a sequence of grid sizes 
Construct and initialize u" 
Define to be the set of equidensity constraints 
Initialize the system of equations V£° 
V-cycle-correction(g;o, u^, C°, A°, V£°) 
Update {x,y) from u*' 
Return (i, y) 



4. Examples of graph drawing layout correction. As previously mentioned, 
the graph drawing problem is of interest for many applications. Therefore, we have 
chosen to demonstrate the abilities of our algorithm for this problem. In this section 
we will present several results of the two-dimensional layout correction algorithm using 
inequality constraints. The set of examples is shown in Figures |4.1| and |4.5| to |4.8[ 
each organized in two columns. The initial and final layouts of the graph are shown 
in the same row, in the left and the right columns, respectively. Note that finalizing 
the "nice graph" representation of these examples is beyond the scope of this work. 
The various "beautifying" procedures used by different applications may, of course, 
be used at the end of our cycles to enhance the visualization results. 

The first example consists of a mesh graph with three holes (Figure [XT] row (a)). 
It is intended to demonstrate that the empty space stays empty and the energy is 



thus kept low. More complicated examples are shown in Figure 4.1 rows (b) and 
(c). The initial optimal positions of the mesh's vertices were randomly changed by 
independent shifts in different directions within a distance d: 



d < 



2hx in example (b) 
Ahx in example (c) 



where hx is the length of a square on the initially taken 32x32 grid, such that the mesh 
size of the graph is actually 2hx. Let us call these meshes Mi and M2, respectively. 
While the correction of AIi looks really nice, two switched vertices at the right-hand 
side of AI2 demonstrate a weak point in our algorithm that certainly must be improved 
by a local "beautifying" procedure, which in general depends on the real application. 
The initial layout (cl) is more complicated than (bl), while the desired final layouts 
are similar. 

A typical example of the energy behavior is presented in Figures |4.2|4.4| These 
figures refer to the mesh example in Figure [XT]- (c). The general energy minimization 
progress is shown in Figure [4^2} In this example the driving routine alternates between 
two grid size V-cycles: each odd V-cycle solves the correction problem for the 16x16 
grid, while even V-cycles improve the previous iterations with the grid 32x32. Figures 



4.3 and 4.4 show the energy behavior of the Window relaxations (without V-cycles) 
for 16x16 grid iterations and alternately 16x16 and 32x32 grid sizes, respectively. 
Clearly, the V-cycle algorithm is more powerful in minimizing the energy than just 
employing the Window relaxations. 



A more complicated example is shown in Figure 4.5 in which the 64x64 mesh 



graph randomly perturbed by vertex shifts (up to 2hx of a 64x64 grid) , compressed at 



the left bottom corner and augmented by 50 randomly chosen edges (Figure 4.5 (a)) 



The final result of the algorithm is presented in Figure 4.5 (b), where all vertices are 



placed almost at their optimal locations (note the different scales of the two figures). 
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We have used 2-FMG cycles with 2 V-cycles at each level as the main driving routine. 
Such a driving routine works with the following grid sizes: 2, 2, 4, 4, 8, 8, 128, 
128, 2, 2, 4, 4, and so forth. After these 2-FMG cycles the total energy was very close 
to its real minimum and additional iterations have only slightly corrected the layout. 
The next experiment consists of the 64x64 compressed mesh with three holes. The 



initial and final layouts are presented in Figures 4.6 (a) and 4.6 (b), respectively 



Two additional examples demonstrate the layout corrections for graphs whose 



vertices have nonequal volumes (see Figures 4.7 and 4.8 1. In both cases the initial 
layout of these graphs was random. 

In spite of the promising results presented above, the algorithm has not yet been 
optimized. However, it is already clear that several parameters must for efficiency be 
kept very small. For example: (1) the number of Window relaxation iterations should 
be fixed between 1 and 3; (2) "optimal enough" in SingleWindowSolver means less 
than 6 iterations and (3) the size of W in SingleWindowSolver is very robust, that 
is, the same results can be obtained with sizes 4x4, 8x8, and 16x16. We have used 
only 4x4 as it runs the fastest. 

5. Conclusions. We have presented a linear time multilevel algorithm for solv- 
ing correction to the nonlinear minimization problem under planar (in)equality con- 
straints. By introducing a sequence of grids over the domain and a new set of global 
displacement variables defined at those grid points, we formulated the minimization 
problem under planar equidensity constraints and solved the resulting system of equa- 
tions by multigrid techniques. This approach enabled fast collective corrections for 
the optimization objective components. We believe that this formulation can open 
a new direction for the development of fast algorithms for efficient space utilization 
goals. Among many possible motivating applications [2J |71 [HI [SI IHl HI ISl IHl EH] we 
focused on the demonstration of the method on the graph visualization problem with 
efficient space utilization demand. 

We recommend this multilevel method as a general practical tool in solving, pos- 
sibly together with other tools, the nonlinear optimization problem under planar 
(in) equality constraints. 
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Fig. 4.2. Energy behavior of the mesh at Figure \4-.l\ (c), when employing complete V-cycles 
with 16 X 16 and 32 X 32 alternately. 
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Fig. 4.3. Energy behavior of Window relaxation iterations (16 x 16 grid) of the mesh at Figure 
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Fig. 4.4. Energy behavior of Window relaxation iterations (16 X 16 and 32 X 32 grids) of the 




Fig. 4.5. Example of the layout of the 64 X 64 mesh with additional random edges (note the 
different scales of the two figures): (a) starting from a compressed and perturbed configuration at 
the bottom-left corner, (b) the resulting picture using V-cycles. 
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Fig. 4.7. Example of the 2D-layout of a graph with nonequal volumes. 
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Fig. 4.8. An example of the 2D-layout of a 5-level binary tree with non-equal vertices. 
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